####################################
# Main effect RDROBUST
####################################

rm(list=ls())

#sink("~/Dropbox/Crime Chile/11_replication/02_rdrobust_crime.txt")

library(Hmisc)
library(ggplot2)
library(stargazer)
library(foreign)
library(rdrobust) 
library(rdd)
library(readstata13)

################
# Prepare data 
################

# read data
d = read.dta13("~/Dropbox/Crime Chile/11_replication/local_crime_data_chile_2022january.dta")   
names(d)

############################
# RDrobust MSE
############################

# NOTE: In tables, we report conventional point estimates and robust confidence intervals and robust p-values.

# Homicide
summary(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster))

pe1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster)$Estimate[1])
pv1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster)$pv[3])
ci1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster)$ci[3,])
oss1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N[1] + rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N[2])
ess1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[1] + rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[2])
bw1 = as.numeric(rdrobust(d$homicide_rate_s,d$margin,all=TRUE,cluster = d$cluster)$bws[1,1])

# Rape
summary(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster))

pe2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster)$Estimate[1])
pv2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster)$pv[3])
ci2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster)$ci[3,])
oss2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N[1] + rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N[2])
ess2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[1] + rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[2])
bw2 = as.numeric(rdrobust(d$rape_rate_s,d$margin,all=TRUE,cluster = d$cluster)$bws[1,1])

# Assault 
summary(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster))

pe3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster)$Estimate[1])
pv3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster)$pv[3])
ci3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster)$ci[3,])
oss3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N[1] + rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N[2])
ess3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[1] + rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[2])
bw3 = as.numeric(rdrobust(d$assault_rate_s,d$margin,all=TRUE,cluster = d$cluster)$bws[1,1])

# Theft
summary(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster))

pe4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster)$Estimate[1])
pv4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster)$pv[3])
ci4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster)$ci[3,])
oss4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N[1] + rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N[2])
ess4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[1] + rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[2])
bw4 = as.numeric(rdrobust(d$theft_rate_s,d$margin,all=TRUE,cluster = d$cluster)$bws[1,1])

# Robbery
summary(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster))

pe5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster)$Estimate[1])
pv5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster)$pv[3])
ci5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster)$ci[3,])
oss5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N[1] + rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N[2])
ess5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[1] + rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[2])
bw5 = as.numeric(rdrobust(d$robbery_rate_s,d$margin,all=TRUE,cluster = d$cluster)$bws[1,1])

# Robbery by surprise
summary(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster))

pe6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster)$Estimate[1])
pv6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster)$pv[3])
ci6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster)$ci[3,])
oss6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N[1] + rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N[2])
ess6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[1] + rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[2])
bw6 = as.numeric(rdrobust(d$robbery_surprise_rate_s,d$margin,all=TRUE,cluster = d$cluster)$bws[1,1])

# TABLE 1
pe = c(pe1,pe2,pe3)
pv = c(pv1,pv2,pv3)
ci_lower = c(ci1[1],ci2[1],ci3[1])
ci_upper = c(ci1[2],ci2[2],ci3[2])
oss = c(oss1,oss2,oss3)
ess = c(ess1,ess2,ess3)
bw = c(bw1,bw2,bw3)
variable = c("Homicide","Rape","Assault")
data = data.frame(variable,pe,pv,ci_lower,ci_upper,oss,ess,bw)

colnames(data) = c("","Point Estimate","Robust P-value","Robust 95% Confidence Interval lower bound","Robust 95% Confidence Interval upper bound","Overall sample size","Effective sample size","MSE bandwidth")
data

stargazer(data, summary=FALSE, rownames=FALSE, out="~/Dropbox/Crime Chile/11_replication/table2.html")

# TABLE 3
pe = c(pe4,pe5,pe6)
pv = c(pv4,pv5,pv6)
ci_lower = c(ci4[1],ci5[1],ci6[1])
ci_upper = c(ci4[2],ci5[2],ci6[2])
oss = c(oss4,oss5,oss6)
ess = c(ess4,ess5,ess6)
bw = c(bw4,bw5,bw6)
variable = c("Theft","Robbery","Robbery by surprise")
data = data.frame(variable,pe,pv,ci_lower,ci_upper,oss,ess,bw)

colnames(data) = c("","Point Estimate","Robust P-value","Robust 95% Confidence Interval lower bound","Robust 95% Confidence Interval upper bound","Overall sample size","Effective sample size","MSE bandwidth")
data

stargazer(data, summary=FALSE, rownames=FALSE, out="~/Dropbox/Crime Chile/11_replication/table3.html")

sink()
